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Abstract. - In the framework of schematic hard spheres lattice models for granular media 
we investigate the phenomenon of the "jamming transition". In particular, using Edwards' 
approach, by analytical calculations at a mean field level, we derive the system phase diagram 
and show that "jamming" corresponds to a phase transition from a "fluid" to a "glassy" phase, 
observed when crystallization is avoided. Interestingly, the nature of such a "glassy" phase 
turns out to be the same found in mean field models for glass formers. 



Gently shaken granular media exhibit a strong form of "jamming" [1-3], i.e., an exceedingly 
slow dynamics, which shows deep connections [4-6] to "freezing" phenomena observed in many 
thermal systems such as glass formers [7] . Although the idea of a unified description of these 
phenomena is emerging [5], the precise nature of jamming in non-thermal systems and the 
origin of its close connections to glassy phenomena in thermal ones arc still open and very 
important issues [8]. 

Here, we discuss these topics in the framework of the Statistical Mechanics of powders 
introduced by Edwards' [9-11] where, to allow theoretical calculations, it is assumed that time 
averages of a system subject to some drive (e.g., "tapping") coincide with suitable ensemble 
averages over its "mechanically stable" states. In particular, we consider a schematic model 
for granular media recently shown [11] to be well described by Edwards' assumption: a system 
of hard spheres under gravity confined on a cubic lattice. In this letter we first show that this 
model subject to a Monte Carlo (MC) "tap dynamics", when crystallization is avoided, has 
a pronounced jamming similar to the one found in experiments [1-3]. We then discuss the 
nature of such a form of jamming by analytically solving Edwards' partition function of the 
system at a mean field level, by use of a Bethe approximation. This approach shows that the 
present model for granular media undergoes a phase transition from a (supercooled) "fluid" 
phase to a "glassy" phase, when its crystallization transition is avoided. The nature of such a 
"glassy" phase results to be the same found in mean field models for glass formers [12, 13]: a 
discontinuous Replica Symmetry Breaking phase preceded by a dynamical freezing point. This 
finding quantitatively confirms early speculations about the structure of jamming in granular 
media [4,14] and clarifies the deep connections with glass formers [5]. 

Our schematic model for granular media is a system of monodisperse hard-spheres (with 
radius — 1) under gravity, constrained to move on the sites of a cubic lattice of spacing oq. 
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Its Hamiltonian is: 



v. = Uhc + mg ^ 



TliZi 



(1) 



where Zi is the height of site i, 5 = 1 is gravity acceleration, m — 1 the grains mass, Ui G {0, 1} 
is an occupancy variable (absence or presence of a grain on site i) and HHciini}) is the hard 
core term preventing the overlapping of nearest neighbors grains [15]. 

N grains are confined in a 3D box of linear horizontal size L and height H (in our MC 
simulations L = 12, H = 2Q, N = 288) between hard walls and periodic boundary conditions 
in the horizontal directions. They are initially prepared in a random loose stable pack and 
are subject to a dynamics made of a sequence of MC "taps" [4]: a single "tap" is a period 
of time, of length tq (the tap duration), where particles can diffuse laterally or upwards with 
a probability pup G [0, 1/2] and downwards with 1 — Pup] when the "tap" is off, grains can 
only move downwards (i.e., pup — 0) and reach a blocked state (i.e., one of their "inherent 
states" [11]) where no one can move downwards without violating the hard core repulsion. The 
parameter p^p has an effect equivalent to keep the system in contact (for a time tq) with a bath 
temperature Tr = mga^i / \a[{l —pup)/Pup] (called the "tap amplitude"). Our measurements 
are performed when the shake is off and the system at rest. Time, t, is measured as the 
number of taps applied to the system. 
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Fig. 1 - Inset The normalized site density correlation, C{t), is plotted, at stationarity, as a function 
of the number of taps, t, for the shown "shaking amplitudes" Tr (tap duration to — lOOMC step per 
site). Main Frame The characteristic relaxation time, r, is shown as a function of Tr. 



The MC tap "dynamics" exhibits large relaxation times which grow as the tap am- 
plitude decreases. In particular, we consider the density correlation function C{t,tw) = 
B{t,ty,)/B{tw,tw), where B{t,t^) = Y.i[{^i(.'t + tw)ni{t^)) - {ni{t + t^)){ni{t^))]. At highTr, 
for tw long enough, C{t,t^) has a time translation invariant behavior, i.e., C(t,tw) = C{t). 
We do not consider here low values of Tr as the model tends to crystallization [16]. The 
relaxation time, T(Tr), plotted in Fig^ has been obtained by fitting C{t) as an exponential 
function at long time t. Even though we are at comparatively high Tr values, an Arrhenius 
law (shown in the picture) or a Vogel-Tamman-Fulcher law fits the data. These properties 
are very similar to those found in more refined models [4, 11, 17] and correspond to recent 
experimental results on granular media [2,3], which are found to exhibit a glassy behavior [7]. 

The nature of the glassy region is, anyway, very difficult to be numerically established 
and further insight can be obtained by analytical treatments. This can be accomplished in 
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Fig. 2 - In our mean field approximation, the hard spheres describing the system grains are located 
on a Bethe lattice, sketched in the figure, where each horizontal layer is a random graph of given 
connectivity. Homologous sites on neighboring layers are also linked and the overall connectivity, c, 
of the vertices isc = A; + l = 5. 

Edwards' approach [9-11] which we now assume to hold in our model, at least as a good 
approximation. This assumption is grounded on some recent results [11] showing that "time 
averages" of the system observables over the "tap" dynamics coincide with those over an 
ensemble average where only "mechanically stable" states, (i.e., those where the system is 
found at rest), are counted. More specifically, the weight of a given state r is [11]: e"''^*^''^ -11^, 
where Tconf = is a thermodynamic parameter, called "configurational temperature", 
characterizing the distribution {Tconf can be related to the shaking amplitude Tr, for instance, 
from the equality between the time average of the energy, e(Tr), and the ensemble average, 
(e)(Tcora/))- The operator 11^ selects mechanically stable states: 11,. = 1 if r is "stable", else 
n,. = 0. The system partition function a la Edwards is thus the following [11]: 

Z^^e-'^^W-n, (2) 

r 

where the sum runs over all microstates r. 

Since the exact calculation of Z for the above lattice model is hardly feasible, we now 
discuss a mean field theory (see [13, 18] and ref.s therein) based on a random graph version of 
such a lattice which is sketched in Fig[21 More specifically we consider a 3D lattice box with 
H horizontal layers (i.e., z £ {1, ...,H}) occupied by hard spheres. Each layer is a random 
graph of given connectivity, fc — 1 (we take k — -i). Each site in layer z is also connected to its 
homologous site in z — 1 and z + 1 (the total connectivity is thus fc + 1). The Hamiltonian is the 
one of eq. plus a chemical potential term to control the overall density. Hard Core repulsion 
prevents two connected sites to be occupied at the same time. In the present lattice model 
we adopt a simple definition of "mechanical stability" : a grain is "stable" if it has a grain 
underneath. For a given grains configuration r = {rii}, the operator 11^ has thus a simple 
expression: 11^ = limx^co Gxp { — ^T~('Edw\ where "H-Edw — ^ni{z),i^ni(z—i),o^ni{z—2),o 

(for 

clarity, we have shown the z dependence in ni{z)). 

The local tree-like structure of our lattice allows to write down iterative equations a la 
Bethe for the probability of fields acting on the lattice sites [19]. In the "cavity method" [18], 
the recurrence equations are found by iteration of the lattice structure where k "branches" 
(i.e., graphs where a root site, denoted by j G {1, ■■■,k}, has only k neighbors) are merged 
to a new site i, leading to a lattice with the same structure as before but with one more 
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site. Three kinds of "branches" exist here: "up" (resp. "down") branches where the root site 
has fc — 1 neighbors on its same layer and one in the upper (resp. lower) layer; and "side" 
branches where the root has k — 2 neighbors on its layer, one in the upper and one in lower 
layer. Correspondingly, three kinds of merging are possible where: an "up" (resp. "down" ) 
branch with root at height z + 1 (resp. z — 1) and fc — 1 "side" branches with root at height 
z merge into a new "up" (resp. "down") branch with root in site i at height z; an "up", a 
"down" and k — 2 "side" branches merge into a new "side" branch. 

The partition function of the new branch ending in site i can be recursively written in 
terms of the partition functions of the merged branches. Define Z^'J'^ and Z^^'^^"^ the partition 
functions of the "side" branch restricted respectively to configurations in which the site i is 
empty or filled by a particle; analogously, z[^'^'' and Zq'^^ (resp. Zq '^^) are the partition 
functions of the "up" branch restricted to configurations in which the site i is filled or empty 

Zq ) are those of the 



with the upper site filled (resp. empty); finally z'^^'^^ and Z^'^"' (resp. 



"down" branch when site i is filled or empty with the upper site empty (resp. filled). For more 
details see [17]. It is convenient to introduce five local "cavity fields" /li*'^'' 



and S^''^"* defined by: 



z 



(t,z) 



= Z 



{^,z) 



0,d 



are more easily written [21]: 



In these new variables the recursion relations 



B(p-mgz) 



fc-2 

11(1 + '^ 



{l + e^a- ) X 



fc-i 



(3) 



fc-i 



f3{p-mgz) -l3h' 



]^(l + e/''^i^-')-i 
i=i 



= (1 + e^^i'" >-''''^' 



From the iterative solution of these equations it is possible to compute the system free energy 
[17, 18]. Figl^l shows the system phase diagram in the plane of the two control parameters, 
{Tconf, Ns), where Ns is the number of grains per unit surface in the box. 

At low Ng or high Tconf a fluid-like phase is found, characterized by a homogeneous Replica 
Symmetric (RS) solution (in Replica Theory terminology) of the recursion equations 0, in 
which only one pure state exists and the local fields are the same for all the sites of the lattice 
(translational invariance). For a given Ns, by lowering Tconf (see Fig.s l3l4|l . a phase transition 
to a crystal phase (an RS solution with no space translation invariance) is found at T,„ [20]. 
Notice that the fluid phase still exists below T,„ as a metastable phase corresponding to a 
supercooled fluid found when crystallization is avoided. 

The above RS fluid solution, however, is not appropriate to describe the high or low 
Tconf region which is dominated by the presence of a large number of local minima of the 
free energy where the fields may fluctuate [18]. This situation is characterized by non-trivial 
probability distributions for the local fields on each layer z: Pz{hu,9u), P^ihs) and Pz{hd,gd). 
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Fig. 3 - The system mean field phase diagram is plotted in the plane of its two control parameters 
Tconf is Edwards' "configurational temperature" and the average number of grains 
per unit surface in the box. At low Ns or high Tconf, the system is found in a fluid phase. The 
fluid forms a crystal below a melting transition line Tm{Ns). When crystallization is avoided, the 
"supercooled" (i.e., metastable) fluid has a thermodynamic phase transition, at a point Tk{Ns), to 
a Replica Symmetry Breaking "glassy" phase with the same structure found in mean field theory of 
glass formers. In between Tm{Ns) and Tk{Ns) a dynamical freezing point, Td{Ns), is located, where 
the system characteristic time scales diverge. 

Within the one-step replica symmetry breaking (IRSB) ansatz of the cavity method [18,22], 
the recursion relations for the fields are replaced by self consistent integral equations for the 
distribution of the local fields: 

X 5[hl-h'i'^^)5{gl~g'i'^^)eM-Pm^K). (4) 

where Ci is a normalization constant, hu'^\ ffu'^' are the local fields defined by the eqs.©, 
A^F!^ is the free energy shifts in the merging process [17,18] and m is the usual IRSB parameter 
to be obtained by maximization of the free energy with respect to it. Analogous equations are 
found for P^{hs) and P^{hd, gd)- We have solved all these equations iteratively, by discretizing 
the probability distributions, until the whole procedure converged. 

A non trivial solution of the IRSB equations appears for the first time at a given tempera- 
ture T]j{Ns), signaling the existence of an exponentially high number of pure states. In mean 
field theory To is interpreted as the location of a purely dynamical transition as in Mode 
Coupling Theory, but in real systems it might correspond just to a crossover in the dynamics 
(see [12,13,23] and ref.s therein). The IRSB solution becomes stable at a lower point Tk, 
where a thermodynamic transition from the supercooled fiuid to a IRSB glassy phase takes 
place (see Fig|21l in a scenario a la Kauzmann with a vanishing complexity of pure states 
(which stays finite ioi Tk <T <Td)- 

The results of these calculations, summarized in the phase diagram of Fig|21 are further 
illustrated in Fig0] in a system with a given number of grains (i.e., a given Ng), the overall 
number density, <&, is plotted as a function of Tconf (here by definition $ = Ns/2{z), where {z) 
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is the average height). The shown curve, ^{Tconf), is the equihbrium function here calculated. 
It has a shape very similar to the one observed in tap experiments [1,3], or in MC simulations 
on the cubic lattice (see also [4]), where the density is plotted as a function of the shaking 
amplitude F (along the so called "reversible branch"). 
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Fig. 4 - For a system with a given number of grains (i.e., a given Ns), the overall number density, 
4> = Ns/2{z} {{z} is the average height), calculated in mean field approximation is plotted as a 
function of Tconf, ^{Tconf) has a shape very similar to the one observed in the "reversible regime" of 
tap experiments and MC simulations of the cubic lattice model for <l>(rr). The location of the glass 
transition, Tk (filled circle), corresponds to a cusp in the function 4>(rco„/). The passage from the 
fluid to supercooled fluid is Tm (filled square). The dynamical crossover point Td is found around the 
flex of ^{Tconf) and well corresponds to the position of a characteristic shaking amplitude F* found in 
experiments and simulations where the "irreversible" and "reversible" regimes approximately meet. 



Summarizing, in the present mean field scenario of a granular medium with particles 
per surface, in general, at high Tconf (i-e. high shaking amplitudes) a fluid phase is located 
(see Figl^l). By lowering Tconf, a phase transition to a crystal phase is found at r,„. How- 
ever, when crystallization is avoided, the fluid phase still exists below as a metastable 
phase corresponding to a supercooled fluid. At a lower point, Tq, an exponentially high num- 
ber of new metastable states appears, interpreted, at a mean field level, as the location of a 
purely dynamical transition, which in real system is thought to correspond just to a dynamical 
crossover. Finally, at a even lower point, Tk, the supercooled fluid has a genuinely thermo- 
dynamics discontinuous phase transition to glassy state. MC simulations of a cubic lattice 
model show indeed the divergence of the relaxation time r at low "shaking" amplitudes (see 
FigQ). The structure of the glass transition of the present model for granular media, obtained 
in the framework of Edwards' theory is the same found in the glass transition of the p-spin 
glass and in other mean field models for glass formers [12, 13]. 
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